Downloading PharmacoSet objects

Let us try to download a PSet onto our local working environment. A table of available PharmacoSet objects can be obtained by using the availablePSets function. Any of the PharmacoSets in the table can then be downloaded by calling downloadPSet, which saves the datasets into a directory of the users choice, and returns the data into the R session. Lets try downloading the GDSC and CCLE PSets.

  availablePSets(saveDir=file.path(".", "PSets"))

  GDSC <- downloadPSet("GDSC", saveDir=file.path(".", "PSets"))
  CCLE <- downloadPSet("CCLE", saveDir=file.path(".", "PSets"))

The basic structure of a PharmacoSet is as follows:

Working with PSets

For most of your work with PSet objects, we recommend you use the accessor functions implemented in the package to access the data stored in the object. For example, here we can use the cellInfo function to pull out the tissue of origin for each cell line in the CCLE dataset. The tissueid column is standard for PSets, and should be present in all objects.

mycol <- c("#8dd3c7","#ffffb3","#bebada","#fb8072","#80b1d3","#fdb462",
           "#b3de69","#fccde5","#d9d9d9","#bc80bd","#ccebc5","#ffed6f",
           "#a6cee3","#1f78b4","#b2df8a","#33a02c","#fb9a99","#e31a1c",
           "#fdbf6f","#ff7f00","#cab2d6","#6a3d9a","#ffff99","#b15928")
pie(table(CCLE@cell[,"tissueid"]), 
    col=mycol, 
    main="Tissue types", 
    radius=1, 
    cex=0.8)

Tissue of origin of cell lines in CCLE study

While we will introduce some of the functions for working with PSets in this tutorial, the documentation of the PharmacoSet object contains descriptions of the available functions to interface with PSet objects:

?`PharmacoSet-class`

PSets can be subsetted by refering directly to the drugs and cell lines you want to keep in the dataset. For example, we can subset the GDSC PSet in the following ways:

print(GDSC)

## Name:  GDSC 
## Date Created:  Wed Dec 30 10:44:21 2015 
## Number of cell lines:  1124 
## Number of drug compounds:  139 
## RNA: 
##  Dim:  11833 789 
## CNV: 
##  Dim:  24960 936 
## Drug pertubation: 
##  Please look at pertNumber(pSet) to determine number of experiments for each drug-cell combination.
## Drug sensitivity: 
##  Number of Experiments:  79903 
##  Please look at sensNumber(pSet) to determine number of experiments for each drug-cell combination.

print(drugNames(GDSC)[1])

## [1] "Erlotinib"

print(subsetTo(GDSC, drugs="Erlotinib"))

## Name:  GDSC 
## Date Created:  Wed Dec 30 10:44:21 2015 
## Number of cell lines:  1124 
## Number of drug compounds:  1 
## RNA: 
##  Dim:  11833 789 
## CNV: 
##  Dim:  24960 936 
## Drug pertubation: 
##  Please look at pertNumber(pSet) to determine number of experiments for each drug-cell combination.
## Drug sensitivity: 
##  Number of Experiments:  323 
##  Please look at sensNumber(pSet) to determine number of experiments for each drug-cell combination.

print(GDSC["YT","Erlotinib"])

## Name:  GDSC 
## Date Created:  Wed Dec 30 10:44:21 2015 
## Number of cell lines:  1 
## Number of drug compounds:  1 
## RNA: 
##  Dim:  11833 0 
## CNV: 
##  Dim:  24960 1 
## Drug pertubation: 
##  Please look at pertNumber(pSet) to determine number of experiments for each drug-cell combination.
## Drug sensitivity: 
##  Number of Experiments:  1 
##  Please look at sensNumber(pSet) to determine number of experiments for each drug-cell combination.

As an exercise, lets manually compute the number of overlaping cell lines and drugs between CCLE and GDSC. Try it yourself!

If you noticed as you completed the previous exercise, cell line and drug names are standardized between PSets downloaded using the downloadPSet function.

PharmacoGx also implements a convenience function to intersect PSets:

common <- intersectPSet(list(CCLE, GDSC), intersectOn=c("drugs", "cell.lines"))
print(common)

Fitting Drug Dose Response Curves

One of the core tasks implemented in PharmacoGx is the fitting of Hill Curve Models to dose-response data. In PharmacoGx, we use the 3 Parameter Hill Slope function as our model of drug response in cancer cell lines: To fit a Hill Slope model to you data, you can use the logLogisticRegression function, as below:

concentrations <- 1/2^seq(0,8) * 1
viabilities <-  c(0, 33.3, 60, 77.8, 88.2, 93.9, 96.9, 98.4, 99.2)

pars <- logLogisticRegression(conc = concentrations, viability = viabilities)

## Warning in sanitizeInput(conc = conc, viability = viability, conc_as_log
## = conc_as_log, : Concentration Values were unsorted. Sorting concentration
## and ordering viability in same order

print(pars)

## $HS
## [1] 1.900066
## 
## $E_inf
## [1] 0
## 
## $EC50
## [1] 0.2813364

Plotting Drug Dose Response Data

Drug-Dose response data included in the PharmacoSet objects can be conviniently plotted using the drugDoseResponseCurve function. Given a list of PharmacoSets, a drug name and a cell name, it will plot the drug dose response curves for the given cell-drug combination in each dataset, allowing direct comparisons of data between datasets.

cells <- c("OCI-AML2","A253","NCI-H1648")
par(mfrow=c(1, 3), pty = "s")
drugDoseResponseCurve(drug="lapatinib", cellline=cells[1], 
                      pSets=list(CCLE, GDSC), plot.type="Fitted", 
                      legends.label="auc_published")
drugDoseResponseCurve(drug="lapatinib", cellline=cells[2], 
                      pSets=list(CCLE, GDSC), plot.type="Fitted", 
                      legends.label="auc_published")
drugDoseResponseCurve(drug="lapatinib", cellline=cells[3], 
                      pSets=list(CCLE, GDSC), plot.type="Fitted", 
                      legends.label="auc_published")

The function drugDoseResponseCurve can also be used to plot your own drug dose response curves, as follows:

concentrations <- 1/2^seq(0,8) * 1
viabilities <-  c(0, 33.3, 60, 77.8, 88.2, 93.9, 96.9, 98.4, 99.2)

drugDoseResponseCurve(concentrations = list("Exp 1" = concentrations), viabilities = list("Exp 1" = viabilities))

## Warning in sanitizeInput(concentrations[[i]], viabilities[[i]], conc_as_log
## = conc_as_log, : Concentration Values were unsorted. Sorting concentration
## and ordering viability in same order

Computing Summary Measures for DDRCs

Often with in vitro Pharmacogenomics data, we want to compare the drug sensitivity of a cell line to some omic feature. For this, we want to summarize the drug dose response curve into a single number representing the sensitivity of the cell line. The IC50 and Area Above the curve are two convenient metrics for quuantifying the observed drug sensitivity. If you noticed above, the drugDoseResponseCurve function computes them by default. In PharmacoGx, they can be computed manually as follows:

concentrations <- rev(1/2^seq(0,8) * 1)
viabilities <-  c(99.2, 98.4, 96.9, 93.9, 88.2, 77.8, 60, 33.3, 0)

print(computeAUC(concentration = concentrations, viability = viabilities))

## [1] 23.68416

print(computeIC50(concentration = concentrations, viability = viabilities))

## [1] 0.2813364

In PharmacoGx, we call these measures sensitivity measures. PSets come with these measures precomputed, and accessible using the sensitivityProfiles function:

head(sensitivityProfiles(CCLE))

Summary Functions

Pharmacogenomics studies often contain many examples of either replicated or missing data. One of the most common tasks in preparing data for statistical or machine learning analysis is aligning your features and labels. To accelerate using Pharmacogenomics data for analysis, PharmacoGx contains two functions which create deduplicated matrices with missing data filled by NAs: summarizeSensitivityProfiles and summarizeMolecularProfiles. They create matrices which are drugs x cell lines and molecular features x cell lines, with any cell lines profiled for only sensitivity or molecular features padded with NA values. For the molecular profiles, the data is returned in a r BiocStyle::Biocpkg("Biobase") ExpressionSet object, while drug sensitivity data is returned as a matrix.

dim(PharmacoGx::summarizeMolecularProfiles(CCLE, "rna"))

## Features  Samples 
##    20049     1061

dim(PharmacoGx::summarizeSensitivityProfiles(CCLE))

## [1]   24 1061

Below, we use the summarizeMolecularProfiles function and ggplot to investigate the distributions of AUC values within CCLE.

library(ggplot2, verbose=FALSE)
library(reshape2, verbose=FALSE)
CCLE.auc <- PharmacoGx::summarizeSensitivityProfiles(CCLE, sensitivity.measure = "auc_recomputed")
melted_data <- melt(CCLE.auc)
NA_rows <- unique(which(is.na(melted_data), arr.ind=T)[,1])
melted_data <- melted_data[-NA_rows,]
p <- ggplot(melted_data, aes(x=Var1,y=value)) +
  geom_boxplot(fill="gray") +
  theme(axis.text.x=element_text(angle=90,hjust=1)) +
  xlab("Drugs") +
  ylab("AAC")
print(p)

Cells response to drugs in CCLE

Consistency of mRNA expression

First we look at the consistency of gene expression profiles across the two datasets. We use the same intersected PSets as above, and calculate Pearson correlations between the expression of common cell lines.

common <- intersectPSet(pSets = list("CCLE"=CCLE, "GDSC"=GDSC), 
                        intersectOn = c("cell.lines", "drugs"), 
                        strictIntersect = TRUE)

## Intersecting large PSets may take a long time ...

common.features <- intersect(rownames(featureInfo(common$GDSC, "rna")),
                             rownames(featureInfo(common$CCLE, "rna")))

GDSC.rna <- PharmacoGx::summarizeMolecularProfiles(common$GDSC, "rna", 
                                       features = common.features)

## Summarizing rna molecular data for:  GDSC

## 
  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |=                                                                |   2%
  |                                                                       
  |===                                                              |   5%
  |                                                                       
  |====                                                             |   7%
  |                                                                       
  |======                                                           |   9%
  |                                                                       
  |=======                                                          |  11%
  |                                                                       
  |=========                                                        |  14%
  |                                                                       
  |==========                                                       |  16%
  |                                                                       
  |============                                                     |  18%
  |                                                                       
  |=============                                                    |  20%
  |                                                                       
  |===============                                                  |  23%
  |                                                                       
  |================                                                 |  25%
  |                                                                       
  |==================                                               |  27%
  |                                                                       
  |===================                                              |  30%
  |                                                                       
  |=====================                                            |  32%
  |                                                                       
  |======================                                           |  34%
  |                                                                       
  |========================                                         |  36%
  |                                                                       
  |=========================                                        |  39%
  |                                                                       
  |===========================                                      |  41%
  |                                                                       
  |============================                                     |  43%
  |                                                                       
  |==============================                                   |  45%
  |                                                                       
  |===============================                                  |  48%
  |                                                                       
  |================================                                 |  50%
  |                                                                       
  |==================================                               |  52%
  |                                                                       
  |===================================                              |  55%
  |                                                                       
  |=====================================                            |  57%
  |                                                                       
  |======================================                           |  59%
  |                                                                       
  |========================================                         |  61%
  |                                                                       
  |=========================================                        |  64%
  |                                                                       
  |===========================================                      |  66%
  |                                                                       
  |============================================                     |  68%
  |                                                                       
  |==============================================                   |  70%
  |                                                                       
  |===============================================                  |  73%
  |                                                                       
  |=================================================                |  75%
  |                                                                       
  |==================================================               |  77%
  |                                                                       
  |====================================================             |  80%
  |                                                                       
  |=====================================================            |  82%
  |                                                                       
  |=======================================================          |  84%
  |                                                                       
  |========================================================         |  86%
  |                                                                       
  |==========================================================       |  89%
  |                                                                       
  |===========================================================      |  91%
  |                                                                       
  |=============================================================    |  93%
  |                                                                       
  |==============================================================   |  95%
  |                                                                       
  |================================================================ |  98%
  |                                                                       
  |=================================================================| 100%

CCLE.rna <- PharmacoGx::summarizeMolecularProfiles(common$CCLE, "rna",
                                       features = common.features)

cors.gexp <- cor(Biobase::exprs(GDSC.rna), Biobase::exprs(CCLE.rna), use="pairwise.complete")
hist(diag(cors.gexp), main="Correlations of gene expression between same cell lines", xlim=c(-1,1))

boxplot(list(Same=diag(cors.gexp), Different = cors.gexp[upper.tri(cors.gexp)]), col=rainbow(5), main="Comparison of gene expression profiles")

Consistency of pharmacological profiles

Before running the correlation analysis, let us examine some examples of drug sensitivity curves which are discordant and concordant across the two datasets, using the functions described above.

common <- intersectPSet(pSets = list("CCLE"=CCLE, "GDSC"=GDSC), 
                        intersectOn = c("cell.lines", "drugs"), 
                        strictIntersect = TRUE)

## Intersecting large PSets may take a long time ...

drugs <- drugNames(common$CCLE)
##Example of concordant and discordant drug curves
cases <- rbind(
  c("CAL-85-1", "17-AAG"),
  c("HT-29", "PLX4720"),
  c("COLO-320-HSR", "AZD6244"),
  c("HT-1080", "PD-0332991"))
par(mfrow=c(2, 2))
for (i in seq_len(nrow(cases))) {
  drugDoseResponseCurve(pSets=common, 
                        drug=cases[i,2], 
                        cellline=cases[i,1], 
                        legends.label="ic50_published", 
                        plot.type="Fitted", 
                        ylim=c(0,130))
}

Consistency of drug response curves across studies

As we can see, there is a range of concordance of drug response between these two datasets. To get a more global picture, we look at the correlations of cell line response between datasets.

GDSC.auc <- PharmacoGx::summarizeSensitivityProfiles(common$GDSC)
CCLE.auc <- PharmacoGx::summarizeSensitivityProfiles(common$CCLE)

corssens <- cor(GDSC.auc, CCLE.auc, use="pairwise.complete")
hist(diag(corssens), main="Correlations of drug sensitivity between same cell lines", breaks=50)

boxplot(list(Same=diag(corssens), Different = corssens[upper.tri(corssens)]), col=rainbow(5), main="Comparison of drug sensitivity profiles")

As we can see, the drug sensitivity data is less concordant between datasets than the gene expression. This presents a challenge for statistical and machine learning, as usually we use the drug response as our target/labels for our algorithms, assuming that they represent a ground truth. As the workshop goes on, we will introduce some suggested techniques to overcome this challenge.

The first method introduces will be the Growth Rate correction techniques implemented in the GRmetrics package.